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Abstract 


Floods are the maximum annual flows of a river and their frequency 
analysis are statistical techniques that allow estimating values 
associated with low probabilities of being exceeded. Such Predictions 
allow the hydrological design and review of hydraulic works. The 
fundamental stage of frequency analysis is the selection and fit, to the 
available data or sample, of a probability distribution function to make 
predictions. The Generalized Extreme Value (GEV) distribution has 
theoretical bases that make ita probabilistic model suitable for floods, 
maximum sea levels, extreme wind speeds, etc. Due to its genesis, the 


GEV processes one data per year, the maximum. This disadvantage 
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has been overcome with the method of the r annual maximum events. 
Such procedure is described and applied to the adjustment of the GEV 
function, through the maximum likelihood method, using the Complex 
algorithm of multiple restricted variables, for numerical maximization. 
Five records of floods with five events per year were integrated from 
hydrometric stations of the Hydrological Region No. 10 (Sinaloa), 
Mexico. Also, two records of maximum sea levels with r = 5 were 
processed, taken from specialized literature, one of them is non- 
stationary since it shows a trend. Results were analyzed and the 
predictions of the GEV adjustment with r annual events, were 
compared against those of the classical L- and LH-moment methods. 
Finally, Conclusions are formulated, which highlight the convenience of 
the described method and point out the advantages of using the 
Complex algorithm as a numerical technique; due to the above, the 
systematic application of the described procedure ¡is recommended. 


Keywords: GEV distribution, maximum likelihood method, Complex 
algorithm, decision and dependent variables, Wald-Wolfowitz test, 


standard error of fit, covariates, predictions. 


Resumen 


Las crecientes son los gastos máximos anuales de un río y su análisis 
de frecuencias son técnicas estadísticas que permiten estimar valores 
asociados con bajas probabilidades de ser excedidos. Tales 
predicciones permiten el diseño y la revisión hidrológica de las obras 
hidráulicas. La etapa fundamental del análisis de frecuencias consiste 
en seleccionar y ajustar, a los datos o muestra disponible, una función 
de distribución de probabilidades para realizar las predicciones. La 


distribución general de valores extremos (GVE) tiene bases teóricas 
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que la convierten en un modelo probabilístico adecuado a crecientes, 
niveles máximos del mar, velocidades extremas de viento, etcétera. 
Debido a su génesis, la GVE procesa un dato por año, el máximo. Esta 
desventaja fue superada con el método de los r eventos máximos 
anuales. Tal procedimiento se describe y aplica al ajuste de la función 
GVE, a través del método de máxima verosimilitud, utilizando al 
algoritmo Complex de múltiples variables restringidas, para la 
maximización numérica. Se integraron cinco registros de crecientes 
con cinco eventos por año, en estaciones hidrométricas de la Región 
Hidrológica No. 10 (Sinaloa), México. Además, se procesaron dos 
registros de niveles máximos del mar con r = 5, tomados de la 
literatura especializada, uno de ellos es no estacionario al presentar 
tendencia. Se analizaron los resultados y se contrastaron las 
predicciones del ajuste de la GVE con r eventos anuales contra las de 
los métodos clásicos de momentos L y LH. Por último, se formulan las 
Conclusiones, las cuales destacan la conveniencia del método descrito 
y señalan las ventajas del uso del algoritmo Complex como técnica 
numérica; debido a lo anterior, se recomienda la aplicación sistemática 
del procedimiento descrito. 


Palabras clave: distribución GVE, método de máxima verosimilitud, 
algoritmo Complex, variables de decisión y dependientes, test de 
Wald-Wolfowitz, error estándar de ajuste, covariables, predicciones. 
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Introduction 


In simple terms, frequency analysis is a statistical technique that allows 
a probabilistic model or probability distribution function (PDF) to be 
fittedtoa sample of n data. Its purpose is to make predictions, beyond 
the data interval used to estimate the fitting parameters of the selected 
PDF. Usually, the intended predictions are greater than the largest of 
the observed events, and therefore frequency analysis consists of the 
study of past events in order to establish probabilities for future 
occurrences. For example, the annual maximum flows of a river are 
studied to obtain the Design Floods, which are maximum flows 
associated with low probability of being exceeded. To summarize, 
frequency analysis methods do not predict future with certainty, but 
generate probabilistic models that explain and make efficient use of 
the extreme events that have occurred in the past (Khaliq, Ouarda, 
Ondo, Gachon 8, Bobée, 2006). 


Frequency analysis consists of four steps: (1) verification of the 
statistical quality of available data or sample. In fact, validity and 
accuracy of the predictions require data to be independent and show 
stationarity, which implies that they have been generated by a stable 
random process over time; (2) selection of a PDF; (3) selection of a 
method for estimating PDF fitting parameters and (4) adoption of 
results or predictions. The standard error of fit is commonly applied, 
which is a quantitative measure between observed and estimated data 
with PDF adjusted (Rao 8 Hamed, 2000; Stedinger, 2017). 
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In relation to the topic of selecting a PDF, the extreme value 
theory isa statistical discipline that develops models that describe the 
unusual nature of data. Being extreme values scarce by nature, their 
predictions involve strong extrapolations, and the extreme value 
theory provides models that allowthem, on an asymptotic basis (Coles, 
2001). 


The Generalized Extreme Value (GEV) distribution has a 
theoretical basis that makes it essential in the analysis of frequencies 
of extreme events, such as river floods, maximum sea levels, extreme 
wind speeds, etc. (Guedes, Soares, €: Scotto, 2004; Khaliq etal., 2006; 
An 8 Pandey, 2007). However, implicit in its own formulation of 
maximum one span or block, which is later discussed, GEV distribution, 
uses only one extreme value per year and therefore two approaches 
have emerged that seek to incorporate more events of each year 
(Coles, 2001; Khaliq et a/., 2006). 


The first one is called peaks over threshold or POT analysis, using 
all events that exceed a threshold value, for example, the lowest value 
of the observed annual maximum. Naturally, the independenae 
condition between the selected events must be met. The second 
approach to using more data is known as the r largest order statistics 
and is about including more information of maximum events in each 
year or block, which meet the independence condition between them. 
In order to fulfill the above, it has been found that r < 5 (Ramesh € 
Davison, 2002). The fundamental advantage of both approaches is to 
use other maximum sample values, at the right end of the GEV 
function, to estimate more accurate or robust predictions, as more data 
are available to quantify their fitting parameters (Tawn, 1988; Coles, 
2001; Khaliq et al., 2006). 
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The objectives of this study focus on the following four: (1) To 
expose the theory behind the application of the GEV distribution, with 
r independent annual maximum events; (2) to describe the maximum 
likelihood method for estimating the three fitting parameters (u1,O,k) of 
the GEV distribution; (3) to integrate five flood records with r = 5, 
within Hydrological Region No. 10 (Sinaloa), Mexico, in order to carry 
out a prediction contrast of GEV distribution, between those obtained 
with the annual r events method and those from the sextile, L moment 
and LH-moment methods. These contrasts are established on the basis 
of the standard error of fit and (4) to process with annual r events 
method, the record of maximum sea levels in Venice, Italy, which is 
not stationary, since it presents a linear upward trend. 


Operational theory 


The GEV distribution 


In frequency analysis of extreme hydrological data, such as floods or 
annual flows, sea levels, wind speeds and maximum rains, extreme 
value distributions play an outstanding role. Its origin is set from a 
sample of random variables (X1, X2, . . . ,Xn), independent and 
identically distributed (i¡d) that have a probability distribution function 
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F(x). In practice, X; are values of the observed process measured 
houny as sea levels, or daily as river flows. 


The statistic defined in Equation (1) is the maximum of the 
process for n units of observation time (block maximum) and its 
probabilistic behavior focuses on development of the extreme value 


theory. 
M,, = max[X,,X>,..., X,,] (1) 


The cumulative probability distribution Mis defined by condition 
¡id and itis (Smith, 1986; Tawn, 1988; Dupuis, 1997; Coles, 2001; An 
8: Pandey, 2007): 


P(M, <x) =P[X,<x%,....,X, <x]1; P(M,_<x) =P(X,<x)*“"P(X, < 
x); P(M, <x) =|[F60]” (2) 


To obtain the behavior of Equation (2) according to n tends to 
infinity and avoid trivial or degenerative solutions, a linear 
normalization of variable M, is performed, based on constants an > O 
and bn; this is: 


q (3) 


Now, extreme value distribution converges to a nontrivial result, 


defined as follows: 
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The asymptotic distribution G(x) must converge in one of the 
three families called Gumbel, Fréchet and Weibull. These three 
distributions are included in the Generalized Extreme Value (GEV), 


whose expression is: 
ae -1/k 
G(x) = exp [1 + a ] (5) 


in which, y, Cd and k are the parameters of location, scale and shape; 
with -o < y < «a, 0 > 0 and -o <k<oo, In addition, it is defined by 
the set Xx: [1+k0=p4)/0]>0). If x> 0 GEV is type Il or Fréchet, with 
no upper limit (u-o/k < x <oo). When x<0 type III or Weibull is 
defined, with upper limit (o < x < u-0/k). Finally, when x =0, in an 
asymptotic sense, the distribution type 1 or Gumbel (-oo < x < 00), 
also called double exponential, is reached, whose equation is: 


cc == (29) 6 


Selecting each year maximums (Equation (1)), the sought 
predictions (Xrr) for an exceedance probability q or a return period Tr 
in years (Tr = 1/q) are obtained based on Equation (5), since GOO) = 
1-q(Coles, 2001): 


Xp = 211 [-m(1- 179 fork + 0 (7) 
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Xy, =4-0*In[—In(1-— q)] fork =0 (8) 


It should be clarifiedthat during the last century GEV distribution 
equation (Smith, 1986; Tawn, 1988; Dupuis, 1997; Hoking € Wallis, 
1997; Rao € Hamed, 2000) was: 


G(x) = exp (- [1 — apo (9) 


For which, the Fréchet model is obtained with x < 0. This change 
of k by —x in Equation (9), has been started, logically by Coles (2001) 
and accepted by Ramesh and Davison (2002), Katz, Parlange and 
Naveau (2002), and by Khaliq et a/. (2006). It has a slight convenience 
of signs when handling the likelihood function and its positivity 
constraints in Equation (13), Equation (14) and Equation (15) below 
and is therefore retained in this study. 


Maximum likelihood method 


According to the referred text (Kite, 1977; Rao 8 Hamed, 2000; Coles, 
2001; Kottegoda 8 Rosso, 2008; Meylan, Favre, € Musy, 2012) the 
maximum likelihood principle is explained in different levels of detail. 


But, in general, it implies the following: Given a sample (X1, X2, . .. 
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¿Xn) of ¡id observations, which follow a distribution Fe with fitting 
parameters 01, 62, . . . ¡0g. Then, by definition, the probability of 
obtaining a value X;, will be: 


P(x¡ <X < x,+ dx,) = fo(x;) : dx; (10) 


with, fe((xi) is the probability density function. As the data are ¡id, the 
probability of obtaining n values of X;, will be the joint probability or 
likelihood function, designated L, and expressed as: 


L(0) = fol)" fol)" fon) = Més fo (5) (11) 


The maximum likelihood method is to find a vector Ó of 
parameters that make maximum L(0) and therefore the probability of 
obtaining the sample (X1, X2, . . . ,Xn). It is often more convenient to 
take logarithms and work with the /logarithmic likelihood function 
(Coles, 2001), ¡.e.: 


(0) =logL(0) = X;=, 108 f4(x;) (12) 


this is acceptable because the logarithmic function is monotonic and 
then the /(0) function reaches ¡ts maximum at the same point as the 
function L(6). 
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Function /(0) of the GEV distribution 


The probability density function of the GEV distribution is as follows 
(Coles, 2001; Ramesh 8 Davison, 2002): 


pen as) 


(0 


where, O > 0, -o < py, k<oo and the x interval is such [1 + k(x-p) / 0] 
> 0, indicated by the + sign outside the square brackets. Any fitting 
parameters combination that violates the previous positivity condition 
implies that at least one of the observed points (x), itis beyond the 
distribution endpoints and then the likelihood function is zero and the 
logarithmiclikelihood function is equal to —oo (Coles, 2001). If you have 
a sample of annual maximum Xi, X2,. . . ,¡Xn, which are independent, 


the logarithmic likelihood function is: 


(0 


lG600 = nin o (142) 82 [14 (202), 82, [1 (222% (14) 


Function /(0) of GEV for r annual events 
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When other annual maximum values are included in the GEV 
distribution fitting, to lookfor more accurate prediction inference; then, 
each year you have r events defined as: x?* > x? >,...,>x” and the joint 
probability density function allows you to define the logarithmic 
likelihood function for n years of rindependent maximum annual values 


each (Coles, 2001; Ramesh € Davison, 2002), which is: l(u,0,k) = —n 


de 
x 4 


r«Ino- (1 + DY, 2-1 In [1 Ae E] 


r 


Ir (E pel 


(0 


The fitting parameters constraints defined in Equation (13) are 
applied to the above expression and [1 +x(x!-— u)/0] must also be 
positive foreach j = 1,..., r. Equation (14) is a special case of Equation 
(15), which is obtained when r = 1. Tawn (1988) presents an equation 
similarto Equation (15), which considers the number of r events to be 
variable each year. 


Function /(0) of GEV for non-stationary series 


Non-stationary processes have characteristics that change 
systematically over time. In climatological and hydrological processes, 


non-stationarity is observed over periods, due to different climatic 
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patterns effects in several months. It is also presentas a trend, mainly 
caused by regional or global climate change (Coles, 2001). 


There is no general extreme value theory for non-stationary 
processes and thus a pragmatic approach is followed to use the GEV 
distribution fitted by maximum likelihood, changing over time or with 
some other covariate, its fitting parameters (Katz eta/., 2002). Forthe 
specificcase of a lineartrend, the location parameter is caused to vary 
over time (Coles, 2001): 


u(t) =P, + Py rt (16) 


The above expression is replaced in Equation (15), in order to 
maximize that function and obtain the four fitting parameters 
((Bo f,,0,x). Other models based on GEV for non-stationary series, with 
the covariate approach, are described by Coles (2001) in his chapter 
6. 


The Complex algorithm 


The maximization of Equation (15) to obtain the optimal fitting 
parameters (u,0a,k) of the sought GEV, must be addressed numerically 
and for this, the Complex algorithm, with multiple variables (2) 
restricted or dimensioned, was selected. Its theoretical approach is as 
follows (Box, 1965): 
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Minimizar F( z,,Z)7,""",Z6) (17) 


Subject to m dependent variables (y), function of decision 
variables (2): 


A EAT E Ta 2) (18) 


Both variables have lower and upper limits of the type <, ¡.e.: 
Zinp £Z¡S Zgup ANA Yins E Y¡ E Ysup- The Complex algorithm is a local 
scanning technique, which is guided solely by what is in ¡its way; its 
background, a brief description of its operational process and its OPTIM 
code, can be found in Campos-Aranda (2003). Another description and 


code of this search method is given in Bunday (1985). 


The primary designations in the OPTIM code are NX and NY, 
which define the number of decision and dependent variables; for the 
case analyzed, three (u, 9, kand n (number of record years), because 
the dependent variables are positivity restrictions, with j = 1,. . ,r 
values per year [1 + k(x! —u)/0]. MI = 500 is the maximum number of 
evaluations of the objective function and NQ = 25 is the number of 
such calculations between printing results. These variables are 
accessed in data reading subroutine. 


An important advantage of the OPTIM code is to allow easy 
access of the limits (L = lower, U = upper), names and initial values of 
the variables, in the above mentioned subroutine, by means of the 
following designations: XL(D, XU(D, XN$S(D, XD, YLO), YUO), YNS(D, 


and Y(J). For the studied case, I varies from 1 to 3, and J from 1 to n 
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(number of record years). Then the convergence criteria FA and FR for 
absolute and relative F deviations are included. The following values 
were used: 0.0002 and 0.00001 respectively. 


The objective function is called F in the OPTIM code and it is 
accessedatthe end ofthe program, logically corresponding to Equation 
(15), with name FO$ =" FLMV”, of the maximum likelihood logarithmic 
function. It is assigned with a negative sign to F, because the Complex 
algorithm minimizes the function (Equation (17)) and it is desired to 
maximize the FLMV, 


Finally, it is clarified that because of the nature of the numerical 
problem raised, we have another variable that is NR = 1,..., 5, which 
is equivalent to each annual sample of levels or flows. Data capture is 
done with two nested cycles, one for í varying from 1 to n or number 
of record years processed and the other for j ranging from 1 to 5, that 
is, r from Equation (15). 


Wald-Wolfowitz test 


This non-parametric test has been used by Bobée € Ashkar (1991), 
Rao € Hamed (2000) and Meylan, Favre 8 Musy (2012) to test 
independence and stationarity in annual maximum flow records (xX;). 
For the studied case, it was proposed to apply the test to the annual 
maximum flow record (r = 1), which should be a sample of random 
values. Wald and Wolfowitz based on the work of Anderson on the 
serial correlation coefficient developed such a test, whose statistic is: 
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R = Ma X¡*Xigq FX" Xg (19) 


When the size (n) of the series or sample (xi) is not small and its 
data are independent, R comes from a Normal distribution with mean 


and variance, given by the following expressions: 


EM =R SS (20) 


n-1 


Var[R] A (21) 
in which: 
S,, = Exp (22) 
Finally, U is calculated, using the equation: 
u= EX (23) 


The U value follows a Normal distribution with a zero mean and 
unit variance, and it can be used to prove the independence of the 
series data with a significance level a, usually 5%. In a two-tailed test, 
the standardized normal variable is Za/2 = 1.96; then, when |U| < 1.96 


the series will be made up of independent values (random sample). 
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Standard error of fit 


It isthe most common indicator(Chai €: Draxler, 2014) for the contrast 
of a probability distribution to real data, established in the mid-1970s 
(Kite, 1977) and it has been applied in Mexico using the empirical 
Weibull formula (Benson, 1962). It will now be applied using the 
Cunnane formula (Equation (25)), which according to Stedinger (2017) 
leads to non-exceedance probabilities (p) approximately unbiased with 
the distributions used in hydrology. The standard error of fit (EEA) 
expression is: 


EEA = Eacaot” 


(n—np) 


(24) 


Xi¡ are the annual maximum levels and flows (r = 1) ordered from low 
to high, whose number is n and X, the estimated maximum flows under 
Equation (7) and the probability evaluated with Equation (25); np = 3 
is the number of GEV fitting parameters. 


(25) 


Search for the optimal and robust solution 
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According to the progress that is made from r = 1 tor = 5, applying 
Equation (15) and finding its maximum with the Complex algorithm, 
more and more data is being processed, first n and finally 5n. This 
modifies the optimal values of the fitting parameters (yu, O, k) and also 
changes the EEA, which is calculated exclusively with the annual series 


of levels or annual maximum flows (r = 1). 


The above, in order to establish an objective comparison, in the 
use of more information of each year processed (robust solution) and 
look for a combination of uy, O and k (optimal solution) leading to an 
EEA, less than or similar to the minimum obtained with the sextile 
adjustment methods (Clarke, 1973; Campos-Aranda, 2001), L 
moments (Hoking € Wallis, 1997; Stedinger, 2017) and LH moments 
(Wang, 1997a,b; Campos-Aranda, 2016). 


Processed data 


North Sea level records 
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Guedes Soares and Scotto (2004) present a graph of five annual 
maximum levels (meters) in the North Sea, in the 1976-1999 period; 
therefore, n = 24. Their approximate values are shown in Table 1. 
These authors, to ensure that they use r independent annual values, 
use the concept of standard storm length or duration (Tawn, 1988), 
defined in 480 hours (20 days), when sampling a record of maximum 
levels taken every three hours. 


Table 1. Five independent annual maximum levels (meters) in the 
North Sea (Guedes Soares €: Scotto, 2004). 


r annual values 


10.38 


10.60 


10.53 10.16 
9.20 8.97 


10.68 9.40 


s sd 
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1989 10.32 10.03 9.80 8.66 
1990 10.46 10.00 8.47 
1991 9.98 9.48 8.98 8.47 8.00 


1992 11.47 9.80 8.46 7.98 7.48 


1994 8.46 7.98 7.20 
1995 11.48 10.98 10.47 9.98 7.98 
1996 10.97 10.00 8.00 700 
1997 12.00 9.33 8.47 7.66 7.48 


1998 8.52 8.40 8.18 8.04 7.88 


Flood records integration 


Five records were processed and integrated based on information 
availableinthe BANDAS system (IMTA, 2002), on CD 1, which is called 
“Maximum Monthly Flows” and includes the day, time and scale reading 
for each of the 12 maximum flows for the analyzed year. The five 
records to be processed from the Hydrological Region No. 10 (Sinaloa), 
Mexico, are presented from greater to smaller basin area. The divide 


and main collector of the hydrometric stations: Huites, Santa Cruz, 
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Jaina, Guamuchil and El Bledal are shown in figure 1 of Campos-Aranda 
(2014), along with 17 other basins in that region. 


For its integration, it was accepted that a period of 15 days, 
assures different meteorological conditions in each flood formation, 
withinthe Hydrological Region No. 10 (Sinaloa), Mexico (Schulz, 1976; 
Campos-Aranda, 2000). 


Extreme event verification was not performed, because it is 
accepted that the hydrometric information contained in the BANDAS 
system (IMTA, 2002) has been checked and purged against field data 


or measurements. 


In order to cause that the records are integrated with independent 
flows, the following process was followed, each year (/): the annual 
maximum flow (xi) is selected and reviewed back and forth of that 
value, if the adjacent monthly flow is at least 15 days apart from the 
date of this first maximum. Monthly flow that does not meet this time 
frame is eliminated. The next flow in magnitude is then sought (xf), 
among the maximum monthly flows available (11, 10 or 9 remaining) 
and the process described is followed. This is the way to integrate the 


five independent annual maximum flows (xj). 
Record of floods in the Huites hydrometric station 


The Huites gauging station in the Fuerte river of the Hydrological 
Region No. 10 (Sinaloa), Mexico, with code number 10037 and basin 
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area of 26057 km?2, started operations in September 1941 and 
concluded in December, 1992 (n = 51), when the construction of the 
Luis Donaldo Colosio Dam began. Withthe process described, its record 
of five independent annual maximum flows, as set out in Table 2, was 


integrated. 


Table 2. Five independent annual maximum flows (m3/s) at the 
Huites hydrometric station in Hydrological Region No. 10 (Sinaloa), 


Mexico. 


rannual values 


1942 2531.0 2 037.6 1 868.8 780.3 427.6 


1943 | 14376.0 | 3283.0 | 1 085.0 416.4 414.2 
1945 1 499.2 1 250.0 268.4 197.4 
1946 1 164.8 445.0 427.5 151.0 
1947 1 127.3 754.8 718.8 634.0 
1948 3 215.0 799.0 623.2 493.5 
1949 | 10000.0 | 2 297.5 895.3 942.4 


1950 ERA O 1 384.0 961.0 439.0 438.9 


1951 677.0 587.5 322.2 37.6 


1952 1 266.0 895.0 355.8 238.3 Md 


1953 1025.0 885.0 68.3 5703 39.0 


1954 954.8 540.1 481.2 406.6 171.4 
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1955 4 780.3 1 069.9 662.0 561.8 


1956 695.7 531.7 493.2 278.6 204.8 
1957 293.0 489.0 380.0 362.0 156.2 
1958 3 010.0 1 045.2 894.0 608.5 307.1 
1959 1 908.0 1 831.0 1 345.5 652.0 544.3 
1961 1 396.3 905.9 831.6 771,2 682.0 
1962 1 620.0 912.0 892.8 501.0 374.5 
1963 2 702.0 1 054.0 980.1 969.2 323.0 
1964 1 319.1 938.5 216.0 145.2 105.0 


1965 1 944.0 1 787.6 491.2 360.0 147.6 


1966 2 420.0 892.2 738.0 688.9 340.0 


1967 2 505.8 1 310.5 1 192.8 462.5 272.6 


1968 1 534.3 1 118.0 1 019.6 703.8 401.6 
1969 1 508.0 736.0 504.0 417.6 300.0 
1970 1 558.0 1 330.0 970.0 328.7 206.0 
1971 2 200.0 1 176.0 592.0 591.2 576.3 
1973 7 960.0 2 256.5 1 546.1 1 380.0 RL 
1974 3 790.0 3 315.0 1 120.0 886.7 587.5 
1975 1 095.0 965.5 532.9 213.2 59.2 
1976 2 677.0 1 350.2 1 211.0 50.4 


1977 1 135.0 622.0 266.9 153.0 120.0 
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1978 4 790.0 1 750.0 1119.2 756.0 475.4 


1979 6 860.0 1 001.0 820.0 480.0 302.0 
1980 1 496.0 1 475.0 1 197.8 660.2 452.0 
1981 4 828.1 2 448.0 2 280.0 2 052.0 1 000.3 
1982 2 450.0 1 571.6 1 085.2 624.0 412.8 
1984 5 580.0 1 623.0 1 132.0 624.5 369.2 
1985 3 585.0 1 250.0 1 121.4 925.0 307.6 


1986 1 348.8 1.329,3 436.8 434.0 193.7 


1987 1 429.2 679.2 310.4 78.7 62.9 


1988 1 866.3 1 494.3 132.6 82.7 47.2 
1990 11 558.6 3 544.2 970.1 815.6 376.9 
1991 2 563.1 2 370.0 1 721.5 1 517.6 1 266.7 


1992 2 025.3 1 564.8 1 348.9 787.3 528.5 


Flood record in the Santa Cruz hydrometric station 


The Santa Cruz gauging station in the San Lorenzo river of the 
Hydrological Region No. 10 (Sinaloa), Mexico, with code number 10040 
and basin area of 8919 km?, began operations in May 1943 and 


concluded its continuous record in December 1980 (n = 38). With the 
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process described, its record of five independent annual maximum 


flows, as set out in Table 3 and Figure 1, was integrated. 


Table 3. Five independent annual maximum flows (m3/s) at the 
Santa Cruz hydrometric station in the Hydrological Region No. 10 
(Sinaloa), Mexico. 


Year rannual values 


IE IG 
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Figure 1. Series of annual maximum flows (r = 1) and stationary 


predictions straight lines at the Santa Cruz hydrometric station in the 
Hydrological Region No. 10 (Sinaloa), Mexico. 


Record of floods in the Jaina hydrometric station 


The Jaina gauging station in the Sinaloa river of the Hydrological 
Region No. 10 (Sinaloa), Mexico, with code number 10036 and basin 
area of 8179 km?, started operations in January 1942 and concluded 
its continuous record in December 1997 (n = 56). With the process 
described, its record of five independent annual maximum flows, as set 


out in Table 4, was integrated. 
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Table 4. Five independent annual maximum flows (m3/s) at the Jaina 
hydrometric station in Hydrological Region No. 10 (Sinaloa), Mexico. 


r annual values 
PO O E 
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Record of floods in the Guamuchil hydrometric 
station 


The Guamuchil gauging stationinthe Mocorito river of the Hydrological 


Region No. 10 (Sinaloa), Mexico, with code number 10031 and basin 
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area of 1645 km, began operations in October 1938 and concluded in 
December 1971 (n = 33), when construction of the Eustaquio Buelna 
Dam began. With the process described, its record of five independent 


annual maximum flows, as set out in Table 5, was integrated. 


Table 5. Five independent annual maximum flows (m3/s) at the 
Guamuchil hydrometric station in Hydrological Region No. 10 


(Sinaloa), Mexico. 


r annual values 


1939 299.0 299.0 164.4 MRE 


1940 254.5 37.2 34.1 
1942 445.0 298.0 161.0 
1943 1 550.0 1 236.4 298.0 


1944 391.8 221 


1945 916.0 336.0 276.0 228.7 


1947 530.0 133.0 15.5 


1948 648.0 548.0 195.9 40.4 2.8 
1949 375.0 145.4 89.2 72.8 34.0 
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1953 1173.0 261.3 213.0 MA 


1954 219.0 115.4 101.2 24.9 
1955 3507.0 189.0 11.4 2.3 
1956 165.0 148.2 76.4 0.4 
1957 526.0 342.0 57.4 8.2 Ma 
AA AA AAA 
1959 1 610.0 374.0 372.8 22.4 
1960 137.0 130.0 94.5 69.3 
1961 524.5 302.0 134.1 28.5 
1962 985.0 524.0 168.0 112.5 


1963 459.5 311,2 21.5 90.7 


1964 390.0 202,1 128,2 MN 


1965 449.0 28Z.> d.L 1.4 


1966 793.9 687.8 3.4 3.e 
1967 TL 03 325.0 105.7 10.6 
1968 200.0 146.2 132.5 24.8 9.3 
1969 312.0 126.6 84.7 30.3 5.3 
- _ 


1971 1 045.0 790.0 175.0 


Record of floods in the El Bledal hydrometric station 
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The El Bledal gauging station in the same name stream of the 
Hydrological Region No. 10 (Sinaloa), Mexico, with code number 10027 
and basin area of 371 km?, began operations in September 1937 and 
concluded in December 1994 (n = 57). With the process described, its 
record of five independent annual maximum flows, as set out in Table 


6, was integrated. 


Table 6. Five independent annual maximum flows (m3/s) at the El 
Bledal hydrometric station in Hydrological Region No. 10 (Sinaloa), 


Mexico. 


r annual values 


1939 eE 204.0 16.5 13.0 
1940 170.0 28.0 19.0 8.4 
1941 118.6 17.1 10.1 4.8 
1942 64.6 33.8 6.4 B.3 


1943 157.2 150.5 32.3 


1944 197.0 1.0 


1945 414.0 253.2 144.0 [AT 
8.4 0.5 0.5 


1947 174.0 


1948 233.0 141.2 106.0 MENE 
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1950 44.1 40.4 22.8 

1951 227.4 68.0 29.6 26.7 

1952 169.0 15.4 15.4 [RE 
Roc Mic cil Mal 

1955 152,0 79.8 319 0 

1956 91.4 31.3 48.1 1.0 

1957 283.0 30.5 28.3 0.3 


1959 132.4 29.3 2.4 £.1 0.3 


ps 


1961 122.0 e PURS 3.8 La 0.5 


1962 1 000.0 108.5 56.5 34.7 
1963 335.0 123.0 61.3 0.4 
1964 258.0 189.8 59.0 12.8 
1965 91.5 56.8 ee pies 
— sd 
1967 325.6 115.4 92.7 46.5 
1968 1 576.0 190.4 70.6 6.1 2.0 


1969 228.0 42.0 24.4 3] 0.3 


1971 276.0 260.0 191.0 
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A 
1974 256.1 98.0 64.2 46.9 25.0 
1975 490.0 73.4 30.6 Ll 0.4 
1976 59.6 25.0 25.0 18.4 6.1 
1978 255.0 21.7 5.6 Ol 01 
1979 283.0 32.3 10.8 0.2 0.1 
1980 226.0 ZL1.L 42.6 1.3 0.2 
1981 210.0 118.0 79.8 1.4 0.3 
1982 67.0 62.5 18.0 11,2 4.5 


1984 160.0 153.1 105.2 


1985 278.0 240.0 95.5 86.5 3.2 


1986 T23:0 20.0 15.1 1355 1,9 
1987 107.2 92.4 11.2 2.4 0.3 
1988 234.3 102.5 80.6 2.8 
_ a 
1990 459.3 332.4 86.2 18.1 
19891 73.8 58.2 35.2 10.3 


1992 181.1 161.9 77.6 12.7 


1994 278.2 143.3 101.0 a 
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Record of levels with trend in Venice 


Table 7 and Figure 2 have the first five annual sea maximum levels, 
out of the 10 that Smith (1986) set out, for Venice, Italy. 


Table 7. Record of five annual maximum levels (meters) for 51 years 
in Venice, Italy (Smith, 1986). 


rannual values Year rannual values 


1931| 103 94 | 1957| 119| 107| 100 
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Figure 2. Series of annual maximum levels (r = 1) and prediction 


straight lines with lineartrend and covariate time (t), in Venice, Italy. 
Results 


Randomization test 


The Wald-W olfowitz test results, defined in equations (19) to (23), for 
each of six processed records (Table 1, Table 2, Table 3, Table 4, Table 
5 and Table 6, columns 2), are as follows: North Sea U = -0.518; 
Huites U = -0.113; Santa Cruz U = -0.827; Jaina U = 1.510; 
Guamuchil U = -1.414; El Bledal U = 0.224, and Venice U = 2.641. 
Therefore, the first six records of annual maximum values (r = 1) are 
random, and that of Venice is not. 


Fittings to annual series (r = 1) 
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Table 8 has concentrated the results of the GEV distribution fittings, to 
each of the annual series of maximum values (r = 1), by means of 
sextile, L-moment and LH-moment methods. These results are the 
fitting parameters (u, O, K) and the standard error of fit (EEA), whose 


minimum value reached in each senes is indicated in parenthesis. 


Table 8. Standard error of fit (EEA) and optimal parameters of the 
annual series (r = 1), in the six records indicated with the three 


mentioned methods. 


¡EE Sea [E 189 a 048 139 -0.3075 


lisa Msacdl 2 Baca 724.345 1 099.256 0.4778 


305.3 814.822 507.620 0.3495 


Jaina 208.5 514.672 390.172 0.4192 


Guamuchil 177.2 373.026 260.903 0.3363 


El Bledal 169.232 118.795 0.3062 


A 3.984 111.438 16.966 -0.1102 


a 


ARE Sea A 184) [A 034 1.187 -0.3040 


Kesede 0 il 541 1178.274 0.4459 
282.2 826.631 402.385 0.4544 
(190.8) 516.752 381.623 0.4353 
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153.4 362.642 250.526 0.3803 


El Bledal 33.4 168.181 123.130 0.2982 


3.636 110.994 16.854 -0.0764 
Record: LH Moments 
IC A A 


North Sea 0.186 10.038 1.106 -0.2391 
(L1) 


Huites (L2) (830.3) 1-509,511 1 598.497 0.3264 


Jaina (L1) 191.6 514.368 387.161 0.4296 


Guamuchil | (152.7) 361.325 249.637 0.3827 
(13) 


El Bledal (32.0) 171.392 116.514 0.3222 
(L4) 


Venice (L1) | (3.587) 111.690 14.241 0.0405 


Santa Cruz | (275.3) 743.090 473.732 0.4044 
(13) 


Guidelines for Complex algorithm application 


At the beginning, we sought to define limits and initial values of fitting 


variables (1, O, K) and dependents or constraints on positivity, which 
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would work for the five series to be processed (r = 1,..., 5) of ea 


record. 


During the first applications of the OPTIM code, several 
difficulties arose, one of which was associated with the limits of the 
variables, which are of the < type. Then for the lower limits you should 
not set zero, but 0.10 forthe location (11) and scale (0) parameters and 
0.01 for the dependent variables. 


With regard to upper limits, it was found advisable to use at least 
twice the calculated values for the GEV parameters u and o, of the 
annual maximum flow series (r = 1), by the method leading to the 
minimum EEA (Table 8). Upper limits, it was never necessaryto change 
them for the OPTIM code to work. A higher value of 100 was adopted 
for dependent variables. 


In relation with the limits of the shape parameter (x), they were 
generally used, as lower -0.50 and higher 0.50; or -1 and 1. Rarely, 
its limit had to be restricted with the aim of obtain results in the order 
of magnitude of the EEA of the first processed record (r = 1). 


Finally, with respect to the initial values of the fitting parameters 
(u, O, K), in general the values near the best GEV fitting worked well 
(Table 8), between sextile, Land LH moments with the series of annual 
maximum flows (r = 1). Table 9 showsthe initial limits and values used 


to process the 30 annual maximum values series. 


Table 9. Limits and initial values adopted forthe fitting parameters 
(y, O, k) in the application of the Complex algorithm. 


Record: 


AU E ONE A 
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North Sea | 0.10 | 100 15 ¡0.10| 100 5 -0.50 | 0.50 | -0.25 


did ds bd 
Gai ho Seis 

ad Boca di hc 

paa pe] 


1000 


0.10 


la le 
El Bledal 0.10 a dei 


The Venice record was processed by replacing Equation (16) in 


0.10 


a ala dd 


Equation (15), now with four fitting parameters. The lower and upper 
limits and the initial value were as follows: f,(0.10,200,100), 
f,(0.10,2,0.75), 0(0.10,30,15) and k(-0.50,0,-0.05). 


When processing some series of certain records, an initial data 
had to be changed so that the Complex algorithm reached a solution 
close to the optimal one and not to a local minimum. For example, in 
the North Sea record, for the series of r = 3, -0.05 was used for the 
initial k. For El Bledal record, the shape parameter value k was limited 
to 0.40. 


Results of the Complex algorithm 


Table 10 shows the main indicators of the Complex algorithm (aC) 
application, such as: The maximum likelihood logarithmic function 
(FLMV), calculated with Equation (15), the standard error of fit (E£EA), 
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evaluated with Equation (23) and the optimum values of the fitting 


parameters (u, O, K). 


Table 10. Results of the Complex algorithm, with each annual series 


of r maximum flows in the indicated records. 


+ FLMV FLMV EA E A 
Initial final [ARA] 
North Sea [E] [pe 9 E] 184 mer 037 1.124 -0.2712 


EEAmín = — 108.8 ——— El 225 Riad 024 0.995 -0.2001 
UnOs E 123.4 cial Lal 0.211 | 10.007 | 1.029 |-0.2164 


158.4 164 0.196 9.993 1.071 -0.2300 


5| 181.3 75.2 198 0.190 10.000 1.118 -0.2753 


Huites 455.1 453.3 5 1363.0 | 1690.25 | 1081.08 0.5713 
EEAmín = 829.5 825.9 a 1186.1 | 1666.73 | 1054.85 0.5490 
830.3 


3| 1185.8 | 1182.5 971.5 | 1759.61 | 1193.36 0.3971 


4| 1515.7 | 1508.9 93 1040.4 | 1781.62 | 1196.44 0.3687 


EEAmín = 2| 552.8 552.1 132 267.3 819.749 | 463.835 0.4040 


5| 1295.4 | 1226.4 | 128 470.6 860.041 | 539.006 0.2255 
Jaina 1| 436.1 435.6 | 127 100.4 505.486 | 372.456 0.4894 
EEAmín = 2| 788.2 784.3 85 230.9 550.876 | 377.100 0.4194 


308 3| 1106.2 | 1098.4 | 120 216.9 550.285 | 376.230 0.4266 


1403.2 | 1394.7 238.3 547.134 | 380.289 0.4131 
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A 1676.0 | 1663.4 134.9 527.104 | 381.679 0.4593 
Guamuchil 243.6 241.9 153.7 | 370.294 | 254.203 | 0.3724 
EEAmín = 448.2 444.9 138.0 | 375.509 | 281.706 | 0.3518 
152.7 


3| 630.5 625.2 103 152.6 338.938 | 267.647 0.3607 


802.1 787.2 119.3 292.374 | 259.547 0.4991 


935.1 917.7 181.7 240.682 | 211.943 0.4994 


El Bledal 374.3 373.5 163.134 | 118.291 0.3673 


EEA nin = 21 684.1 | 676.5 | 110 | 33.8 | 144.617 | 116.827 | 0.3473 
32.0 31 9514 | 931.5 | 167 | 42.0 | 131.231 | 108.428 | 0.3998 


4| 1202.2 | 1159.1 | 167 71.2 115.991 99.014 0.4000 


1367.1 | 1326.9 140.1 93.270 76.361 0.3999 


It is also cited in column 1 of Table 10, the minimum EEA 


obtained (Table 8) with the annual maximum series (r = 1), with one 
of the three methods applied: sextile, L moments and LH moments. 
Against this value, the values obtained by the Complex algorithm, 
shown in column 6 (robust solution), must be compared. The adopted 
EEA values are shown shaded in Table 10, are the minimum values 


achieved by the Complex algorithm. 


The results for the five series of the Venice, Italy record have 
been concentrated in Table 11; in which it is observed that the best 
fitting was achieved with r= 2 and an EEA value close to the minimum 
of the LH-moment method of 3.587 meters (Table 8). 


Table 11. Results of the Complex algorithm, with each annual series 
of r maximum levels in the Venice, Italy record, considering linear 


trend. 
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Initial final 


MEAN 


E 


Bs 


Discussion 


EEA contrasts 


Based on the results in Table 10, the following is deducted: (1) for the 
North Sea record, that the maximum likelihood method with r = 1, 
equals the minimum EEA, achieved with the L—-moment method. (2) in 
the records of the Huites and Santa Cruz stations, the best option is 
obtained with r = 5, in the first one, and r = 2, in the second one. In 
Santa Cruz, an EEA less than the minimum (L3) of the LH-moment 
method is obtained. (3) the following are defined as best options in the 
records of the stations Jaina, Guamuchil and El Bledal: r= 1,r=4 and 


r= 1, leadingto an E£EA less than the minimum obtained with the L- 
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moment method in Jaina, and the LH-moment method in the other two 
stations. 


Predictions of the Complex algorithm 


Selecting, in each of the six records processed, the r series that led to 
the smaller EEAs, the predictions obtained and set out in Table 12 are 
contrasted against those of the method that led to the smaller EEAs in 
Table 8, in order to establish which ones will be adopted, because they 
are more severe or critical; they appear shaded. 


Table 12. Predictions for the indicated return periods, with each 
annual series of r independent maximum flows in the indicated 
records. 

Record: Met. Return periods (Tr), in years 
r 
North Sea | ml 11.46 11.97 12.46 | 12.75 | 12.97 | 13.35 | 13.46 


(meters) 1 11.42 11.93 12.44 | 12.74 | 12.99 | 13.41 | 13.54 


11.31 11.83 12.38 | 12.72 | 13.02 | 13.56 | 13.75 
11.33 11.84 12.38 | 12.72 | 13.01 | 13.52 | 13.70 


4 11.35 11.87 12.42 | 12.75 | 13.03 | 13.53 | 13.70 


11.37 11.88 12.35 | 12.67 | 12.92 | 13.33 | 13.46 
4603 6821 10523 | 14113 | 18593 | 33829 | 43286 
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4256 6642 11562 | 17380 | 25997 | 65651 | 97674 


4123 6355 10869 | 16112 | 23758 | 57976 | 84967 
4206 6099 9457 | 12906 | 17428 | 34194 | 45434 


4 4178 5976 9089 | 12213 | 16228 | 30607 | 39953 
5 4432 6600 10581 | 14805 | 20501 58214 
Bac ad isis od Baal Nica bl bad 
Ha 1770 2450 3624 | 4796 | 6297 | 11637 | 15087 
2 1776 2521 3851 5225 | 7035 | 13801 | 18370 
3 1772 2428 3499 | 4516 | 5762 9870 | 12349 
4 1804 2500 3644 | 4738 | 6086 13308 


5 1822 2440 3387 4232 5214 8174 9817 


mL 1324 1975 3168 4431 6133 | 12745 | 17364 


1330 2034 3386 | 4882 | 6975 | 15670 | 22108 
1338 1962 3091 | 4271 5843 | 11832 | 15945 


1341 1972 3120 4328 5945 | 12162 | 16465 


1337 1959 3077 | 4240 | 5783 | 11615 | 15593 
Ap” 1351 2032 3307 | 4684 | 6569 | 14116 | 19526 


Guamuchil 1252 1928 2613 | 3503 6744 | 8883 
(m8/s) s88l 1266 1934 | 2607 | 3474 6592 | 8628 
932 1342 2042 2735 | 3615 6702 | 8672 


872 1268 1949 2628 | 3497 6576 | 8560 


A a e er 
AAA Era ld id Bss Loaf li 
AAA 
II 


1112 1470 | 2719 3512 


Tecnología y ciencias del agua, ISSN 2007-2422, 12(4), 163-218. DOI: 10.24850/j-tyca-2021-04-04 


2021, Instituto Mexicano de Tecnología del Agua 
Open Access bajo la licencia CC BY-NC-SA 4.0 
Tecnología y 


Ciencias¿Agua 


(https://creativecommons.org/licenses/by-nc-sa/4.0/) 


It is clarified that the predictions shown for the five series in 
Venice, Italy, correspond to the end of the historical period, ¡.e., using 
the covaniate t with a value of 51, in Equation (16), when applied in 
Equation (7). 


It is clarified that due to the extent shown (33 < n <57) by the 
five records of floods processed for the Hydrological Region No. 10 
(Sinaloa), Mexico, their reliable predictions may span up to the 100- 
year retumn period (Tr); however, those of Tr = 500 and 1 000 years 
are presented, to observe the degree of dispersion that such robust 


predictions show at both extreme recurrence intervals. 


For the Santa Cruz record, in return periods (Tr) greater than 50 
years, slightly larger predictions are shown in Table 12 with the LH 
moments method, than with robust solutions, of the annual maximum 
r events method. The same happens in the Venice record, but only in 
the return periods (Tr) of 500 and 1000 years. 


Moreover, for the records of Huites (r = 5), Santa Cruz (r = 2) 
and Guamuchil (r = 4), more severe or critical predictions were 


obtained with the fitting ofthe GEV and the method of r annual events, 


Tecnología y ciencias del agua, ISSN 2007-2422, 12(4), 163-218. DOI: 10.24850/j-tyca-2021-04-04 


2021, Instituto Mexicano de Tecnología del Agua 
Open Access bajo la licencia CC BY-NC-SA 4.0 
Tecnología y 


Cienciasz Agua (https://creativecommons.org/licenses/by-nc-sa/4.0/) 
ls! 


than those achieved with classical methods using a maximum annual 


value. The above is remarkable in the Guamuchil record. 


Finally, in the remaining records (North Sea, Jaina and El Bledal) 
with r = 1, the maximum likelinhood method equals or reduces the EEA, 
achieved by the classical methods of L and LH moments. At the Jaina 
station, itis feasible to adopt the results of the series with r = 5, if a 
robust solution is desired, with a slightly larger EEA. 


Conclusions 


Fitting the GEV distribution, using the maximum likelihood method and 
taking r annual maximum events, is a procedure that allows to use 
more information on the observed extreme values, and therefore 
makes more reliable predictions, as its fitting parameters are being 
estimated more accurately. To summarize, itleads to a robust solution, 


from a statistical point of view. 


The seven numerical applications described showed that 
maximizing the logarithmic likelihood function, using the Complex 
algorithm, did not present computational difficulties. 1t also finds 
optimal solutions that improve the achieved fitting with the classic L- 
and LH-moment methods, which use only the annual maximum value. 
This was the case, in four of the five flood records processed (Table 
10). 
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Based on the latest numerical application, relative to the five 
maximum sea levels in Venice, Italy, the simplicity of processing non— 
stationary records, which present a linear trend, is observed. This is 
due to the ease of increasing decision variables, changing their limits 
and/or initial values and modifying the objective function in the 
Complex algorithm. 


With regard to predictions, which varied from a return period of 
5 to one of 1000 years, solutions achieved with the Complex algorithm 
generally reported larger or more critical values. This was remarkable 
inthe records of Guamuchil, Jaina and El Bledal. The exception was the 
record ofthe Santa Cruz station, inthe high periods of return, for which 


slightly lower predictions were obtained (Table 12). 


Taking into account that it does not involve great difficulty to 
integrate records of independent floods with five maximum values per 
year, and that the implementation of the Complex algorithm is simple, 
it is recommended to systematically apply the fitting of the GEV 
distribution, with rannual events in the analysis of floods frequencies. 
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